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ABSTRACT 


The  computer  generation  of  pseudo-random  variates  from 
the  beta  distribution  is  discussed.  Presented  are  three  exact 
methods  applicable  for  parameter  values  p > 1 and  q > 1. 

All  three  methods  use  rejection  from  regions  defined  by  the 
location  of  the  points  of  inflexion.  The  methods  differ  in  the 
amounts  of  set-up  required,  with  higher  set-up  times  generally 
resulting  in  lower  marginal  generation  times.  Generation  times 
are  compared  to  three  methods  in  current  use:  those  of  J'dhnk, 

Fishman,  and  Ahrens  and  Dieter.  Over  a wi'J~  range  of  parameter 
values  marginal  generation  times  are  re  *bout  50%. 

Substantial  set-up  time  savings  and  moderate  marginal  time 
savings  are  obtained  compared  to  Forsythe's  approach  as  recently 
developed  by  Atkinson  and  Pearce. 

KEY  WORDS 

Beta  Distribution 
Simulation 

Random  Variate  Generation 
Rejection  Methods 
Monte  Carlo 
Distribution  Sampling 
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1.  INTRODUCTION 


t 


The  pseudo-random  generation  of  values  from  the  beta  distri- 
bution is  considered  in  this  paper.  The  standardized  beta  distribu- 
tion has  density  function 

fg(x)  - xP  1 (l-x)q  1/S(p,q)  0 < x < 1 


where  the  beta  function  is  defined 

S(p,q)  - f1  uP_1  (l-u)'1'1  du 
0 


Since  the  more  general  form  of  beta  variates  over  the  range 
[a,b]  may  be  obtained  by  the  simple  transformation  y a + x(b-a) , 


only  the  standard  form  is  considered  here.  This  paper  considers 

parameter  values  p > 1 and  q > 1.  For  these  parameter  values 

lim  fQ(x)  * 0,  lim  fQ(x)  * 0,  and  there  is  a unique  mode  at 
x-0  3 x-1  6 


-i  x - (p-1)  / (p+q-2) . 


f 


1.1  Application 

Over  the  interval  [a,b]  the  beta  distribution  is  used  as  a 
model  to  describe  random  phenomenon.  A common  example  is  to  repre- 
sent activity  times  In  PERT  networks  (See,  e.g.,  Clark  (1962)). 

The  beta  distribution  is  well-suited  to  this  purpose,  as  it  can 
assume  a wide  variety  of  shapes.  Schmeiser  and  Deutsch  (1977) 
Include  beta  generators  in  their  discussion  of  general  process 
generators,  due  to  the  distribution's  ability  to  assume  a wide 
range  of  skewness  and  kurtosls  values.  In  particular,  the  Bernoulli 
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trial  and  uniform  distributions  are  special  cases  and  the  gamma 
and  normal  distributions  are  limiting  cases. 

The  beta  distribution  also  arises  in  a number  of  theoretical 
ways.  The  kC^  order  statistic  of  a sample  of  size  n of  uniform 

(0,1)  values  has  density  function  fQ(x;  k,  n-k+1) . Ramberg  and 

p 

Tadikamalla  (1977)  use  this  result  to  generate  order  statistics 

from  distributions  having  a closed-formed  inverse  cdf.  Any  method 

of  generating  from  the  beta  distribution  with  fractional  p and  q 

is  also  a method  for  generating  from  other  related  distributions. 

In  particular,  if  X is  distributed  fQ0r,  -7) , then  Y =*  (w/v)X/(l-X) 

p l Z 

has  an  F distribution  with  v and  w degrees  of  freedom  (Patel,  Kapadia 

1/2 

and  Owen  [1976]).  If  X is  beta  (1/2, w/2),  then  Y =*  s[wX/(l-X)]  has 
a Student's  t distribution  with  w degrees  of  freedom,  where  s is  a 
random  sign.  The  Pearson  Type  VI  distribution 

f(y)  * yw/((i-+y)v  S(wH,  v-w-i))  y > 0 

may  be  generated  using  Y » X/(l-X)  where  X has  a beta  (w+1,  v-w-1) 
distribution.  These  and  other  applications  are  discussed  in  Johnson 
and  Kotz  (1970) . 

1.2  Survey  of  the  Literature 

Two  approaches  for  generating  pseudo-random  beta  variates 
have  appeared  in  the  literature:  those  based  on  specific  proper- 
ties of  the  beta  distribution  and  those  based  on  rejection  from 
the  density  function.  They  are  discussed  in  turn. 

Three  methods  based  on  specific  properties  of  the  beta  distri- 


bution have  been  proposed  and  used.  When  p and  q are  integer  the 
first  method,  due  to  Fox  (1963),  returns  the  pth  order  statistic 
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of  a sample  of  size  q+p-1  from  the  uniform  distribution  as  a beta 
(p,q)  observation.  Rather  than  sorting  to  determine  the  order 
statistics,  the  order  statistics  may  be  generated  more  efficiently 
by  using  the  technique  of  Schucany  (1972)  if  p is  close  to  q+p-1  or 
the  technique  of  Lurie  and  Hartley  (1972)  if  p is  close  to  one.  The 
second  method,  due  to  JVhnk,  is  based  on  X * Y^P/(Y^P  + Z^'*) 
being  beta  (p,q)  if  Y^P  + Z^"^  < 1,  where  Y and  Z are  independent 
uniform  (0,1)  random  variables.  Both  Johnk's  method  and  the  third 
method  are  valid  for  all  p > 0 and  q > 0.  The  third  method  is 
based  on  the  result  that  X * Y / (Y  + Z)  has  a beta  distribution  if 
Y and  Z are  independent  gamma  random  variables  with  shape  parameters 
ex  - p and  * q,  respectively.  While  any  exact  method  for  gener- 
ating the  gamma  variates  could  be  used,  the  method  presented  by 
Fishman  (1973)  is  considered  here:  the  product  of  a U-shaped  beta 

variate  (p,  q < 1)  and  an  exponential  variate  are  added  to  an  Erlang 
variate  to  obtain  the  appropriate  gamma  variate.  (Other  techniques 
for  generating  gamma  random  variates  may  be  found  in  McGrath  and 
Irving  (1973),  Ahrens  dnd  Dieter  (1974),  Wallace  (1974),  Whittaker 
(1974),  and  Fishman  (1976)). 

Three  other  methods  for  beta  generation  based  on  rejection  via 
comparison  with  the  density  function  have  been  proposed.  The  first 
method,  applicable  when  p > 1 and  q > 1,  is  the  general  rejection 
technique  consisting  of  the  generation  of  uniformly  distributed 
points  (x,v)  in  the  rectangle  formed  by  (0,0),  (0,  f(p-l,  p+q-2)), 
(1,  f(p-l,  p+q-2)),  and  (1,0)  and  the  acceptance  of  x as  the  beta 
variate  if  v < f(x).  (See,  e.g.,  Lewis  (1975)).  A better  fitting 
region  than  the  rectangle  was  proposed  by  Ahrens  and  Dieter  (1974) . 
Their  algorithm  BN  generates  uniformly  distributed  points  (x,v)  in 
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a region  based  on  the  normal  distribution.  For  the  parameter  values 
p > 1 and  q > 1,  the  probability  of  acceptance  is  increased  considerably 
compared  to  the  first  rejection  method.  Recently  Atkinson 
and  Pearce  (1976)  proposed  a method,  but  did  not  specify  an  algorithm, 
based  on  the  work  of  Forsythe  (1972).  An  elaborate  set-up  step,  in- 
volving numerical  integration,  was  used  to  allow  efficient  generation 
in  terms  of  marginal  generation  time  per  variate.  A fourth 
method,  efficient  when  at  least  one  parameter  is  less  than 
one,  is  the  "switching"  algorithm  given  by  Atkinson  and  Whit- 
taker (1976).  This  rejection  method,  which  is  based  on  decom- 
posing Che  density  function  into  a product  of  two  functions, 
is  not  efficient  when  both  parameter  values  are  greater  than 
one. 

1.3  Motivation  and  Direction  of  the  Research 

This  paper  is  motivated  by  the  need  for  techniques  capable  of 
generating  pseudo-random  beta  variates  faster  than  existing  techniques 
allow.  The  commonly  encountered  case  of  p > 1 and  q > 1 is  considered. 
In  many  situations  requiring  computer  generated  beta  variates,  many 
observations  are  obtained  for  fixed  values  of  p and  q.  It  is  in 
these  situations  that  the  three  techniques  presented  here  will  be 
of  greatest  benefit.  By  spending  some  additional  time  in  initiali- 
zation, the  marginal  time  per  variate  is  reduced. 

Two  tactics  are  employed  to  reduce  the  number  of  time  consuming 
operations  necessary  to  generate  each  variate.  The  composition- 
rejection  method  (see,  e.g. , Jansson  (1966)),  using  only  triangular 
and  rectangular  regions  from  which  it  is  easy  to  generate  random 
points,  is  combined  with  a preliminary  rejection  and/or  acceptance 
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step  which  reduces  the  number  of  times  the  beta  density  function 
must  be  evaluated.  Variations  of  these  tactics,  combined  with  some 
standard  practices,  produce  three  techniques  which  have  various 
initialization,  marginal  time,  and  memory  requirements.  The  three 
methods  are  developed  in  section  2.  In  section  3 they  are  compared 
to  existing  methods. 


> 


\ 
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2.  DEVELOPMENT  OF  THE  TECHNIQUES 

Developed  in  this  section  are  three  related  techniques  for 
generating  beta  variates  for  p > 1 and  q > 1.  All  three  techniques 
generate  a point  (x,v)  distributed  uniformly  over  a region  defined 
by  a function  t(x)  which  envelopes  f(x)  (i.e.,  t(x)  > f(x)  for  all 
xe[0,l]).  The  value  x is  delivered  as  the  next  observation  if 
v < f(x).  Tocher  (1963)  proves  that  such  an  approach  is  valid. 

This  rejection  approach  is  modified  in  all  three  techniques 
presented  here.  A function  b(x)  is  defined  such  that  b(x)  < f(x) 
for  xe[0,1].  Then  x is  delivered  as  the  next  observation  if 
or  if  v < f(x).  By  defining  b(x)  as  a piece-wise  linear 
--on,  substantial  computation  is  saved  by  the  capability  of 
accepting  x based  on  comparisons  with  b(x)  rather  than  f(x),  since 
f(x)  involves  time  consuming  operations.  Note  that  the  use  of 
b(x)  reduces  the  computation  required  to  accept  or  reject  x,  but 
does  not  change  the  probability  of  acceptance  or  rejection. 

It  should  be  noted  that  rejection  techniques  do  not  need  to 
consider  the  normalizing  constant;  in  this  case  S(p,q)  However, 
3imply  ignoring  S(p,q)  causes  scaling  problems.  Therefore,  following 
Ahrens  and  Dieter  (1974),  define  P - p-1,  Q - q-1,  R - P + Q,  and 

f(x)  - (x/P)P((l-x)/Q)QRR  0 < x < 1 (1) 

which  is  proportional  to  fg(x),  does  not  involve  8(p,q)  and  has 
no  scaling  problems  since  evaluation  at  the  mode  always  yields 
a value  of  one. 


All  three  techniques  presented  in  this  paper  use  the  result 
that  the  points  of  inflexion  of  the  beta  distribution  lie  at 
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[P  ± (pq/(r-i))1/2]/r 


if  the  values  lie  between  zero  and  one  and  are  real.  The  density 
function  is  concave  between  the  points  of  inflexion  and  convex 
at  all  other  points  in  [0,1].  In  addition  for  p > 1 and  q > 1, 
the  mode  is  located  at  P/R  which  lies  between  the  points  of 
inflexion.  (See,  e.g.,  Johnson  and  Kotz  (1970)). 

For  any  P > 0 and  Q > 0,  define  the  constants 


x 


1 


X 


2 


x2(l-x2)/ (P-Rx2) 
(PQ/(R-l))1/2)/R 


if  x2  - 0 
otherwise 

if  real  and  in  [0,1] 

otherwise 


x3  - P/R 


1 1I 

l(P  + (PQ/(R-1))A,*)/R 


+ x4(1-x4)/(P-Rx4) 


(2) 


if  real  and  in  [0,1] 
otherwise 

if  x.  ■ 1 
4 

otherwise 


The  mode  is  represented  by  x^.  If  there  are  inflexion  points, 
they  are  represented  by  x?  and  x4*  If  there  are  no  points  of 
inflexion,  then  x2  ■ 0 and  x4  ■ 1.  x^  and  x^  are  the  points 
at  which  the  lines  tangent  to  the  density  function  at  x0  and 
x,  cross  the  X axis.  If  x»  “0,  then  x.  * 0 and  if  x.  “ 1 , 
then  Xj  • 1.  Figure  A illustrates  the  location  of  these  points. 


J 


Figure  A About  Here 


2.1  Modified  Ahrens  and  Dieter  Technique 

The  Ahrens  and  Dieter  (1974)  algorithm  BN  is  repeated  here  for 
convenience . 


1.  Set  P-ep-l,  Q«-q-l,  R P + Q,  L + R In  R, 

U - P/R,  a - 0.5/R1/2 

2.  Take  an  observation  s from  the  standard  normal  distribu- 
tion and  form  x *■  scr  + y. 

3.  If  x < 0 or  x > 1 go  to  step  2. 

4.  Generate  u ^ U(0,1). 

5.  If  In  u > P In  (x/P)  + Q In  ((l-x)/Q)  + L + .5s2  go  to  step  2. 
Otherwise,  deliver  x as  the  beta  (p,q)  observation. 

2 1/2 

Here  the  enveloping  function  is  t(x)  * exp(-s  /2)  where  s * 2(x-P/R)R 
The  modification  proposed  is  to  sometimes  bypass  the  time  consuming 
step  5 through  the  use  of  the  preliminary  acceptance  function 


b^  (x) 


I (x-x1)/(x3-x1) 

(X5~X)/ (Xj-X^) 


0 < x < x^ 

X^  < X < x^ 

x^  < x < x,. 

x5  < X < 1 


as  illustrated  in  Figure  B.  Proposition  1 forms  the  basis  for  the 
modified  Ahrens  and  Dieter  technique  BNM  proposed  here. 


Figure  B About  Here 
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Proposition  1:  For  x^,  and  x^  defined  as  in  equation  (2), 

b^(x)  < f(x)  for  all  xe[0,l]. 

Proof:  The  proof  is  trivial  for  0 < x < x^  and  x_  < x < 1,  since 
b^(x)  ■ 0 < f(x)  over  these  intervals.  Now  consider  the  interval 
[x^f  x^].  Let  v(x)  ■ f(x7)  + f ' (x.,)  (x-x2)  , the  tangent  to 
f(x)  at  the  inflexion  point  x^.  Define  i(x ) as 


g(x) 


y(x),  xetx.^,  x2] 
f(x),  xe(x2,  x^]. 


Then  g(x)  is  differentiable  everywhere  on  [x^,  x^]  since 


lim 

h-*0 


g(x2+h)-g(x2) 


g(x-+h)-g(x.) 

lia  r — - f ' UJ  • 

h-0"  h 2 


Also  g'(x)  is  monotone  non-increasing  on  [x^,  x^] , first  on  [x^,  x2 
by  linearity  of  y(x),  and  second  on  (x2>  x^]  by  concavity  of  f(x) 
(Roberts  and  Varberg  (1973)).  Hence  g(x)  is  concave  on  [x^,  x^] . 
Moreover,  y(x)  < f(x)  on  [x^,  x2]  since  f(x)  is  convex  in  this 
region  and  y(x)  is  the  line  of  support  at  x?.  It  follows  that 
g(x)  < f(x)  for  all  xe[x^,  x^].  Now  for  x£[x^,  x^],  let 
x * Ax^  + (l-X)x^  for  0 < A < 1.  Then 


f(x)  > g (x)  > Ag(x^)  + (1-A)  g(x3) 

■ Ab^x^  + (1-X)  b1(x3) 

- b^(x)  by  linearity  of  b^*)- 

Similar  logic  is  applicable  to  the  interval  [x3>  x^],  completing 
the  proof  of  Proposition  1. 

The  Ahrens  and  Dieter  technique  may  now  be  modified  as  follows 

a.  Add  to  step  1 the  calculation  of  x^  and  x,.  (which  require  x2 

and  x,  as  intermediate  values)  and  note  that  x.  ■ u- 
4 3 
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b.  Between  steps  4 and  5 add  step  4b  as  follows: 
„ 2 

If  u exp  (-s  /2)  < b^(x) , deliver  x. 
Otherwise  go  to  step  5. 

No  other  changes  are  necessary. 


2.2  The  "Two-Point"  Technique 

A second  method  for  generating  beta  variates  is  developed 
in  this  section.  The  normal  enveloping  function  of  the  Ahren's 
and  Dieter  approach  is  replaced  in  this  section  with  a piece-wise 
linear  function  t^(x) . Calculations  necessary  to  determine  t^(x) 
allow  the  use  of  a better  fitting  preliminary  acceptance  function  b^(x) . 
functions  t^(x)  and  b^(x)  require  f(x)  to  be  evaluated  at  x^  and 
x^  in  step  1.  This  method  will  be  referred  to  as  the  two-point 
technique  2P.  It  is  illustrated  in  Figure  C. 


Figure  C About  Here 


Consider  the  enveloping  function  t^(x)  defined  as 


r 


ci(x)  " { 


x f(x2)/x2 
1 
1 


i^(l-x)  f(x4)/(l-x4) 


° < X < x2 

X,  < X < x3 

x3  < x f x4 

X.  < X < 1 

4 - 


Proposition  2:  For  x^,  x^,  and  x^  defined  x in  equations  (2), 

t^(x)  > f(x)  for  all  xe[0fl]. 

Proof:  Consider  0 < x < x2>  Since  f(x)  is  convex  over  the  interval, 

fUl-XHOJ+tocj)  < (l-X)f(O)  + Xf^) 


The 


1 


-4 


11 


f(Xx^)  < AfCx^) 

Since  x ■ Xx£  in  this  interval, 
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to  f(x)  at  x?,  hence  b^(x)  < f(x)  by  convexity  of  f(x).  For  xefx^.x^], 
b,,(x)  is  a cord  Co  f(x)  and  b2(x)  < f(x)  by  concavity  of  f(x)  in  this 
interval.  It  follows  that  b9(x)  < f(x)  for  all  xetx^.x^].  Similar 
logic  applies  to  prove  that  b2(x)  < f(x)  for  all  xetx^.x^],  thus 
proving  the  proposition. 

Just  as  the  function  b^(x)  can  be  used  to  accept  x without 
comparison  to  f (x) , the  normal  density  function  can  be  used  to 

2 

reject  x without  comparison  to  f(x).  Define  $(x)  ■ exp(-2R(x-x2)  )• 
The  point  (x,v)  may  be  rejected  if  v > <p(x)  since  <p(x)  > f(x)  for 
all  x£ [0,1] , as  proved  by  Ahrens  and  Dieter  (1974). 

The  two-point  technique  proceeds  as  follows: 

1.  Calculate  x^,  x^,  x^,  and  x^  using  equations  (2) 

and  f(x0)  and  f(x^)  using  equation  (1). 

2.  Generate  a value  x with  probabilities  proportional 
to  t^ (x) • 

3.  Generate  a value  v distributed  U(0,  C^(x)). 

4.  If  v < b^(x)  deliver  x. 

5.  If  v > <p(x)  go  to  step  2. 

6.  If  v < f(x)  deliver  x,  otherwise  go  to  step  2. 

Actual  implementation  of  the  algorithm  proceeds  best  by 
considering  each  interval  separately,  due  to  the  fragmented 
definitions  of  t^(x)  and  b2(x). 


2.3  The  "Four  Point"  Technique 

The  third  technique  for  beta  variate  generation  is  given 
here.  Based  on  evaluating  f(x)  at  x^,  x2>  x^,  and  Xj,  the  "four 
point"  technique  4P  uses  the  better  fitting  enveloping  function 
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f 


xfCx^/x^ 

f(x1)-Kx-x1)  [f(x2)-f(X;L)  I /(x2-  Xx) 


1 


f(x5)-Kx5-x)  [f(x4)-f(x5)J  /(x5-x4) 
^ (1-x) f (x5)/(l-x5) 


0 < x < x. 


Xx  < X < x2 
x2  < X < x3 


X < X < X. 
j - 4 


X.  < X < xc 
4 - 5 

x5  c X < 1 


as  illustrated  in  Figure  D.  Proposition  4 states  that  t2(x)  is  indeed 
an  enveloping  function. 


Figure  D About  Here 
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Proposition  4:  For  x^,  x2>  x^,  x4,  and  x^  defined  in  equations  (2), 

t2(x)  > f(x)  for  all  xe[0,l]. 

Proof:  The  proof  follows  the  logic  of  the  proof  of  Proposition  2, 

but  considers  the  intervals  [O.x^],  [x^x^ , [x4,x,.]  and  [x.,1] 
individually.  The  logic  is  not  repeated  here. 

For  the  implementation  of  the  4P  algorithm  t^(x)  involves 
trapezoidal  regions  for  the  intervals  [x^,x2]  and  [x4,x^].  Values 
of  x can  be  generated  from  each  trapezoid  through  generating  from 
a rectangle  and  a triangle  with  probabilities  proportional  to 
their  areas.  However,  since  each  rectangle  is  entirely  below 
f(x)  (see  fig.  D)  all  the  points  (x,v)  generated  with  the  rectangle 
probabilities  will  be  accepted  and  no  rejection  test  need  be  per- 
formed in  this  case.  Likewise,  for  all  the  points  (x,v)  generated 
with  probabilities  proportional  to  the  inner  rectangle  areas  for 
the  [x2>x2J  and  [x2>x4]  regions  no  rejection  test  is  needed.  Thus 
th*  implementation  of  the  four  point  technique,  also  given  in  the 
Appendix,  differs  from  the  two  point  technique  in  that  four  regions 
of  zero  rejection  probability  are  created. 
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3.  COMPUTATIONAL  EXPERIENCE 

In  this  section  the  three  methods  developed  in  section  2; 
designated  8NM,  2P  and  4P;  are  compared  to  the  Ahrens  and  Dieter 
method  BN,  the  ratio  of  gamma  random  variates  method  RG,  and 
Jdhnk's  method  JK.  Section  3.1  discusses  the  implementation 
of  the  comparison,  and  section  3.2  discusses  the  results  of 
the  comparisons . 

3.1  Implementation 

All  timing  comparisons  were  made  on  the  SMU  CDC  CYBER  72 
using  the  KRONOS  2.1.2-A  operating  system.  Each  of  the  six 
algorithms  were  implemented  in  FORTRAN  using  the  FTN  compiler 
at  optimization  level  0.  The  U(0,1)  values  were  obtained  from 
RANF,  the  relatively  fast  generator  built  into  FTN. 

Each  method  was  coded  as  efficiently  as  is  reasonable.  The 
Ahrens  and  Dieter  methods,  both  original  and  modified,  were 
implemented  using  the  Kinderman  and  Ramage  (1976)  algorithm 
for  generating  normal  deviates,  as  this  appears  to  be  the 
fastest  method  available.  (The  initial  use  of  another  generator 
resulted  in  two  to  three  times  inflation  of  computation  time.) 

The  ratio  method  RG  was  implemented  using  the  approach  described 
in  Fishman  (1973).  The  source  listings  for  all  programs  used 
to  obtain  the  computational  results  are  available  from  the  authors. 

For  each  combination  of  p,  q and  method,  four  samples  of 


size  1000  were  generated.  Never  is  the  ratio  of  standard  devia- 
tion of  average  time  to  average  time  greater  than  .01.  The  least 
significant  digit  shown  is  correct  to  within  one  or  two  units. 
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3.2  Computational  Time  Comparisons 

Symmetric  distributions  lead  to  the  fastest  generation  times 
for  all  six  methods.  Table  1 shows  marginal  generation  times  for 
the  six  methods  over  a wide  range  of  parameter  values  p = q,  cor- 
responding to  the  symmetric  case.  Method  JK  deteriorates  rapidly 
as  p ■ q increases,  due  to  the  increasing  probability  of  rejection. 
Method  RG  is  quite  competitive  for  small  integer  values,  but  requires 
two  beta  generations  using  JK  for  fractional  p ■ q and  requires 
p + q uniform  observations  for  all  values  of  p and  q,  making  RG 
inefficient  for  large  and/or  fractional  parameter  values.  Methods 
BN  and  BNM  perform  best  for  large  parameter  values,  due  to  the 
asymptotic  normality  of  the  beta  distribution. 

It  is  clear  from  Table  1 that  the  three  methods  presented  in 
this  paper  reduce  marginal  computation  time  substantially.  Method 
BNM  reduces  the  time  of  BN  by  282  for  p ■ q > 100,  with  less  reduc- 
tion for  smaller  parameter  values.  For  p and  q less  than  two,  the 
4P  reduces  to  the  logic  2P  and  hence  the  times  are  very  similar. 

For  larger  values  of  p and  q,  4P  is  substantially  faster. 

Table  1 can  be  used  to  compare  the  three  methods  of  this  paper 
with  Forsythe's  method.  As  applied  to  the  beta  distribution  by 
Atkinson  and  Pearce  (1976),  Forsythe's  method  requires  632  to  852 
of  the  marginal  time  needed  by  the  BN  method  for  parameter  values 
ranging  from  two  to  seven.  Thus,  over  this  range,  Forsythe's 
method  appears  to  be  about  as  fast  as  BNM,  and  somewhat  slower 
than  2P  and  4P.  However,  the  Forsythe  method  requires  a set-up 
step  equal  to  about  1000  marginal  generation  times.  This  is 
considerably  more  than  required  by  any  of  the  six  methods  compared 


in  Table  1. 
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1.  Marginal  Generation  Times  (in  y Seconds) 
for  Symmetric  Beta  Distributions 


METHOD 


p-q 

JK 

RG 

BN 

BNM 

2P 

4P 

1.01 

1.2 

2.0 

1.8 

1.8 

.37 

.39 

1.2 

1.4 

2.2 

.85 

.80 

.39 

.41 

1.5 

1.9 

2.3 

.77 

.69 

.42 

.43 

2 

3.4 

.56 

.77 

.64 

.45 

.46 

3 

11.5 

.62 

. 74 

.60 

.31 

.27 

5 

* 

.77 

.72 

.55 

.35 

.26 

10 

it 

1.1 

.69 

.51 

.44 

.28 

100 

* 

7.1 

.67 

.48 

1.2 

. 45 

1000 

* 

★ 

.67 

.50 

3.6 

.97 

* Not  calculated 

due  to  excessive  time  requirement. 

A 
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The  set-up  times  for  the  six  methods  vary  widely.  Both  JK 
and  RG  require  no  set-up  time,  since  neither  makes  use  of  constants 
based  on  p or  q.  The  BN  method  requires  calculation  of  a natural 
logarithm  and  a square  root.  The  modified  version,  BNM,  requires 
an  additional  square  root.  The  2P  and  4P  methods  require  varying 
amounts  of  set-up  time  ranging  from  no  higher  order  operations 
to  the  evaluation  of  f(x)  two  and  four  times,  respectively. 

Table  2 shows  generation  times  for  the  six  methods  when  the 
set-up  logic  is  performed  for  each  observation  generated.  Methods 
JK  and  RG  are  unchanged,  since  they  require  no  set-up.  For 
parameter  values  less  than  two,  2P  and  4P  are  fastest  for  both 
set-up  and  marginal  generation  times.  For  parameter  values  greater 
than  two,  the  largest  change  from  Table  1 is  for  4P,  with  2P,  BNM 
and  BN  requiring  decreasing  set-up  time.  For  these  cases 
symmetry  was  not  used  as  a computational  advantage.  Therefore  the 
times  shown  are  also  indicative  of  the  results  for  asymmetric 
distributions. 

It  can  be  seen  from  Tables  1 and  2 that  the  additional  set-up 
time  of  4P  usually  results  in  lower  marginal  times  than  BN  or  BNM. 

The  number  of  values  which  must  be  generated  to  break-even  ranges 

from  less  than  one  for  small  values  of  p and  q to  about  five  for 
p ■ q - 10  and  to  about  60  for  p »■  q - 100.  Thus  4P  is  fastest 

except  for  large  values  of  p and  q combined  with  small  sample  sizes. 

Table  3 shows  marginal  execution  times  for  all  six  methods 
for  a range  of  parameter  values  corresponding  to  asymmetric  dis- 
tributions. Several  points  can  be  made  concerning  Table  3.  First 
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2.  Generation  Times  (in  u Seconds) 
Including  One  Time  Set-Up  For 
Symmetric  Beta  Distributions 


METHOD 


p-q 

JK 

RG 

BN 

BNM 

2P 

4P 

1.01 

1.2 

2.0 

2.2 

2.0 

.69 

.83 

1.2 

1.4 

2.2 

1.1 

.99 

.69 

.84 

1.5 

1.9 

2.3 

1.0 

.86 

.68 

.81 

2 

3.4 

. 56 

1.0  ’ 

.96 

.87 

1.0 

3 

11.5 

.62 

1.0 

1.0 

1.7 

2.7 

5 

* 

.77 

.99 

.98 

1.7 

2.6 

10 

* 

1.1 

.99 

.95 

1.8 

2.7 

100 

* 

7.1 

.98 

.91 

2.6 

3.1 

1000 

* 

* 

.99 

.94 

5.1 

3.6 

* Not  calculated  due  to  excessive  time  requirement. 


■A 
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3.  Marginal  Generation  Times  (in  u Seconds) 
for  Asymmetric  Beta  Distributions 


Method 


V 

1 

JK 

RG 

BN 

BNM 

2 P 

4P 

1.01 

1.5 

1.4 

2.2 

.91 

.82 

.44 

.44 

2 

1.7 

1.3 

.93 

. 74 

.48 

.47 

5 

3.4 

1.4 

1.2 

1.1 

.91 

.44 

10 

6.0 

1.5 

1.6 

1.5 

1.8 

.77 

100 

* _ 

4.5 

4.6 

4.7 

16.1 

5.7 

1.5 

2 

2.5 

1.5 

.78 

.69 

.44 

.46 

5 

6.4 

1.5 

.82 

.70 

.47 

.33 

10 

* 

1.7 

.97 

.87 

.81 

.42 

100 

* 

4.7 

2.5 

2.6 

6.7 

2.3 

2 

5 

12.0 

.67 

.78 

.63 

.43 

.35 

10 

* 

.83 

.89 

.76 

.68 

.40 

100 

* 

3.8 

2.0 

2.0 

4.9 

1.6 

5 

10 

* 

.93 

.74 

.58 

.45 

.28 

100 

* 

3.9 

1.4 

1.3 

2.7 

.84 

10 

100 

it 

4.0 

1.1 

1.0 

1.8 

.62 

* Not  calculated  due  to  excessive  time  requirement. 
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note  that  RG  is  fast  for  small  integer  parameters  and  is  fairly  slow 
for  fractional  parameter  values.  Second,  note  that  BNM  saves  10-20% 
compared  to  BN,  except  for  large  parameter  values.  Finally,  note 
that  4P  dominates  the  other  methods  for  all  parameter  values  shown, 
except  for  p ■ 1.01  and  q =*  100  where  RG,  BN,  and  BNM  are  all 
faster. 

3.3  Other  Comparisons 

Only  the  speeds  of  the  six  methods  were  considered  in  section 
3.2.  Some  comments  concerning  other  criteria  are  given  here.  Since 
all  are  exact  techniques  no  discussion  of  accuracy  is  given. 

The  storage  requirements  of  the  method  varies  widely.  In 
order  of  ascending  storage  requirements  the  methods  rank  JK, 

RG,  2P,  4P,  BN  and  3NM.  This  ranking  includes  the  use  of  the 
Kinderman-Ramage  normal  generator.  At  the  expense  of  computa- 
tional speed,  the  use  of  another  normal  generator  requiring  less 
core  could  place  BN  and  BNM  before  2P  and  4P. 

Generality  is  the  other  important  criterion  for  selection 
of  generators.  Fox's  method,  using  order  statistics,  is  valid 
only  for  integer  parameter  values.  Only  JK  and  RG  are  theoretically 
valid  for  any  p > 0 and  q > 0.  The  other  methods  apply  only  for 
p > 1 and  q > 1. 


•4 
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4.  SUMMARY  AND  CONCLUSIONS 

The  three  rejection  methods  BNM,  2P  and  4P  for  generating 
beta  variates  for  p > 1 and  q > 1 have  been  developed  and  compared 
to  three  existing  methods  JK,  RG,  and  BN.  As  a group,  the  three 
new  methods  dominate  the  three  existing  methods  in  terms  of 
generation  time.  In  particular,  the  best  overall  performer  of 
the  existing  methods,  BN,  is  dominated  by  BNM  for  all  parameter 
values  in  terms  of  both  marginal  time  and  set-up  time.  However, 

BNM  does  not  dominate  either  2P  or  4P.  Except  for  some  insigni- 
ficant cases,  method  4P  dominates  method  2P,  due  to  (•)  fitting 
the  tails  better  than  t^(*)* 

In  conclusion,  method  4P  is  fastest  except  in  the  following 
three  cases: 

1)  If  parameter  values  are  integer,  p + q < 20,  and  the 
parameter  values  change  for  each  observation  generated, 
method  RG  is  fastest.  Note  that  in  this  case  Fox's  order 
statistics  method  is  appropriate. 

2)  If  the  standard  deviation  is  less  than  approximately 
.02  BNM  is  fastest,  due  to  the  good  fit  of  the  normal 
density  function  and  the  corresponding  heavy  tails. 

(A  method  6P  would  be  effective  here.) 

3)  If  p and  q are  greater  than  two  and  the  parameter  values 
change  for  each  observation  generated,  then  BNM  is  fastest 
due  to  its  relatively  fast  set-up  compared  to  2P  and  4P. 


Methods  2P  and  4P  may  be  modified  in  a variety  of  ways.  2P  could 
include  two  additional  regions  having  zero  probability  of  rejection. 
Both  2P  and  4P  could  use  various  other  t(*)  and  b(0  functions.  In 


22 


particular,  b2(x)  could  be  made  to  better  fit  f(x)  in  4P  by  making  use 
of  f'Cx^)  and  f'(x_).  Finally,  methods  corresponding  to  6P,  8P,  ... 
could  easily  be  developed  based  on  the  results  of  section  2.  While 
these  modification  would  be  very  useful  for  large  parameter  values,  the 
added  complication  would  probably  not  be  worthwhile  for  most  users. 
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